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suppression in relativistic heavy ion collisions. 
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1. Introduction 

The theory of strong interactions describes the fundamental force responsible for 
the stability of nuclei and the particles that themselves constitute nuclear matter, 
e.g. the proton, neutron, etc. The theory describes the interaction and properties 
of the most elementary constituents of matter currently known: quarks and gluons. 
The quantum theory which describes the interaction of color-charge carrying quarks 
and gluons is called Quantum Chromodynamics (QCD). QCD predicts that when 
the temperature of nuclear matter is increased above a certain threshold, strongly 
interacting matter undergoes a phase transition to a "new" state of matter referred 
to as the quark-gluon plasma (QGP). Through this phase transition the degrees of 
freedom change from color-neutral hadrons to color-charged partons which are no 
longer confined to exist only inside color-neutral hadrons. We note that calling this 
phase of matter "new" QGP is somewhat of a misnomer, since it is thought to be 
one of the oldest forms of matter in the universe, presumably present a millionth 
of a second after the birth of the universe in the Big Bang. QCD predicts that in 
the high temperature deconfined QGP phase, the interactions between quarks is 
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color-screened much like electric charge is screened in an electromagnetic plasma. 
In addition, the QGP is also characterized by the restoration of a fundamental QCD 
symmetry, which is broken at low temperatures, called chiral symmetry. 

Due to non-linearities present in QCD there is as of yet no comprehensive ana- 
lytic understanding of the theory. For difficult non-perturbative questions, of which 
there are many in QCD, practitioners instead utilize numerical calculations of the 
QCD partition function on discrete lattices called lattice QCD. As a result there 
are several supercomputer centers around the world utilizing state-of-the-art com- 
putational resources in order to aid our understanding of fundamental questions 
in QCD. As mentioned, lattice QCD is a quite reliable tool particularly for ad- 
dressing the non-perturbative aspects of QCD. The phase transition itself is highly 
non-perturbative and, in the immediate vicinity of the transition, one cannot rely 
on weak-coupling perturbation theory techniques written solely in terms of either 
hadronic or quark and gluon degrees of freedom. Many modern finite tempera- 
ture lattice QCD results have achieved a high-degree of precision and the results 
obtained are consistent between different groups using a variety of lattice discretiza- 
tion methods. Unfortunately, the determination of a precise transition temperature 
is impossible since, as we will discuss below, the QCD phase transition is a smooth 
crossover and not a sharp first or second order phase transition. 

One can, however, see deconfinement of the partonic degrees of freedom through 
a rapid rise in thermodynamic quantities, e.g. energy density and pressure, around 
T ~ 160 MeV, which indicates the liberation of quark and gluon degrees of free- 
dom. At temperatures of approximately 300 MeV and above, the entropy density 
found in lattice studies has an approximately 20% deviation from the entropy den- 
sity of a non-interacting gas of quarks and gluons. This deviation is well-described 
by resummed weak-coupling QCD calculations, indicating that at these tempera- 
tures partonic degrees of freedom in the form of quark and gluon quasiparticles are 
the relevant degrees of freedom. At such temperatures it is clear that color-electric 
screening becomes relevant. The transition from hadronic degrees of freedom to 
quark and gluon degrees of freedom is a smooth process occurring between 160 and 
300 MeV and, therefore, as mentioned previously, there is no well-defined transition 
temperature. There is, however, a chiral symmetry restoration transition tempera- 
ture that can be clearly identified. Lattice studies find the chiral symmetry restora- 
tion temperature to be at T ~ 154(9) MeV. In the rest of this work, and in the 
majority of the literature, the chiral symmetry restoration temperature is referred 
to as the QCD critical temperature, T c . 

Since the 1970's quarkonia, which are bound states of heavy flavor quarks (charm 
or bottom) and their antiquarks, have had an essential role in testing QCD, espe- 
cially at zero temperature. Since the seminal work of Tetsuo Matsui and Helmut Satz 
in 1986p^ the behavior of quarkonia in high-temperature QCD has been considered 
a tell-tale signal for deconfinement and QGP production in heavy ion collisions. The 
original idea behind this proposal was that color screening in a deconfined plasma 
would reduce the binding of quarks, and thus prohibit the formation of heavy quark 
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bound states. This "melting" of quarkonia in the QGP would manifest itself in 
the suppression of the quarkonia yields relative to the production without a QGP 
formation. Quarkonium melting at high temperature has been often quoted to be 
the QGP thermometer (see Fig. Q] for a cartoon rendition of the idea), serving as 
an essential tool to measure the plasma temperature. Various states with different 
masses have different sizes and thus are expected to be screened and dissociate in 
the QGP at different temperatures. 

As we discuss in detail in this review, the ex- 
istence and the dissociation of quarkonia in the 
hot QGP medium is far from being this sim- 
ple. The three leading complications are that 
(i) there are more than just color screening ef- 
fects affecting the binding of a particular state 
in the QGP (there is for instance Landau damp- 
ing), (ii) there exists a poorly understood non- 
perturbative intermediate temperature region 
just above the chiral critical transition (the 
aforementioned domain 160 — 300 MeV) where 
identifying the relevant degrees of freedom, and 
thus their relevance for quarkonia must be care- 
fully considered, and (hi) there exist purely nu- 
clear initial state effects, collectively called cold 
nuclear matter effects (e.g. nuclear absorption, 
parton energy loss, and modifications of the par- 
ton distribution function), which must also be 
accounted for. 

During the last two decades, quarkonia in a finite temperature medium have 
been studied using lattice QCD, but the results, more precisely the interpretation 
of the results, have not been without controversy. Effects of the hot medium are 
imprinted on the spectral function, the quantity which contains all the informa- 
tion in a given heavy quark-antiquark channel. Therefore, in principle, determining 
spectral functions at different temperatures should give us insight into the effects 
of a QGP on the bound states. For instance, the disappearance of a well-defined 
peak into the background of continuum states would be a clear indication for the 
melting of a state. This allows us to define what is referred to as a dissociation tem- 
perature (as we discuss below, we would like to draw attention to inconsistencies in 
the literature in how dissociation temperatures are precisely defined) . In this review 
we discuss state-of-the-art lattice QCD results relevant for in-medium quarkonium 
properties, including screening and time- and space-like heavy meson correlation 
functions. Both of these are directly connected to the spectral functions that we 
seek and enable a determination of the static quark-antiquark interaction potential. 
We also briefly review what we have learned about quarkonia from these results. 

The formation, behavior, and modification of the properties of the quarkonia has 



Fig. 1. A cartoon of the QGP ther- 
mometer provided by the sequential 
melting of quarkonia taken from He I. '1 
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also been addressed in phenomenological models, commonly referred to as potential 
models. We would like to highlight though, that the semantics is somewhat mis- 
leading, because potential models need not rely on phenomenological constructions 
of the heavy quark interaction potential, as was done in early studies. The potential 
can be derived directly from the QCD Lagrangian. At zero temperature, this was 
accomplished some time ago already. In recent years, finite temperature field the- 
ory derivations of the real part and, the previously ignored, imaginary part of the 
heavy quark potential have been completed. We will refer to these QCD-derived ap- 
proaches as potential models even, and discuss them in more detail below. There are 
currently limitations to what questions such potential results can answer, making 
other phenomenological models necessary. Therefore, when educated guesses for the 
potential (based on high temperature QCD, and/or lattice QCD results for the free 
energy) are used to analyze quarkonia at finite temperature, or dynamical modeling 
is utilized to bridge from static spectral functions to their dynamical evolution in a 
cooling and expanding medium, we refer to these as phenomenological models. 

In addition to the theoretical work, the QGP is studied in heavy-ion collision 
experiments at different collision energies at the SPS (CERN), RHIC (BNL) and 
the LHC (CERN). The QGP has been produced and found to behave as a nearly 
perfect fluid, although the detailed properties and behavior of this matter are still 
under scrutiny. These experimental programs dedicate significant effort to measur- 
ing quarkonium production. An important experimental observable is the nuclear 
modification factor, Raa, which is the observed particle yield in nucleus- nucleus 
collisions (AA) relative to scaled proton-proton collisions (pp) where QGP forma- 
tion is not expected. Quarkonium suppression has been observed in all of the above 
experiments. The ground state charmonium J/ip suppression was a key signature 
of deconfinement already at the SPS. More recent experiments have led to the 
understanding that not all of the suppression is caused by the QGP: quarkonium 
suppression observed in proton-nucleus collisions (pA) indicates suppression due to 
cold nuclear matter effects. Understanding the wealth of experimental results re- 
quires understanding the medium modification of quarkonia and disentangling hot- 
and cold-matter effects. In view of the wealth of data available from RHIC and LHC 
on charmonium and bottomonium suppression, with particular emphasis on their 
ground states, the J /if) and the T, we present our current theoretical understanding 
of quarkonia. 

We note in closing that recently the question of quarkonia at finite tempera- 
ture has been addressed using string theory-based AdS/CFT methods. Due to the 
nascent nature of these studies, we feel it is still premature to draw firm conclusions 
from these studies concerning our understanding QCD and, in particular, heavy ion 
collision data. Since space for this review is limited and these studies fall outside 
the immediate area of expertise of all three authors, we don't discuss them further 
in this review; however, we are encouraged that there may be an additional non- 
perturbative method for addressing the behavior of quarkonia in the quark gluon 
plasma. 
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Our review has the following structure: In Section 2 we present the state-of- 
the-art in theoretical methods for addressing questions concerning quarkonia in the 
QGP. In particular, we discuss relevant results from lattice QCD and from the effec- 
tive field theory approach to potential models. In Section 3 we address quarkonium 
phenomenology, including quarkonia in a viscous QGP, the status of dynamical 
models, sequential suppression, and cold nuclear matter effects. In Section 4 we 
present our conclusions and outlook for future studies. 

2. Theoretical Methods 

2.1. QCD at Non-zero Temperature 

The properties of hot strongly-interacting matter can be studied using first princi- 
ples field theoretical methods. The QCD partition function can be written as a path 
integral over the gluon and quark fields living in Euclidean space-time with a finite 
time direction r < 1/T, with T being the temperature P The effects of the medium 
are encoded in the boundary conditions on the gluon and quark fields: gluon fields 
obey periodic boundary conditions, while quark fields obey anti-periodic boundary 
conditions. This formalism is commonly known as the imaginary time formalism. 
Dynamical properties of hot strongly interacting matter can be obtained from time- 
dependent correlation functions. In the imaginary time formalism one can calculate 
correlation functions along the imaginary time axis. These can be related to real 
time correlation functions by analytic continuation. However, in many cases it is 
simpler to formulate the problem directly in real time. This can be done using the 
real time formulation of thermal QCD, in which the time runs in the complex plane, 
first along the real axis then coming back from infinity to — ijT ? 

The imaginary time formalism allows one to calculate the path integral numeri- 
cally using lattice discretization (see Refs.|4]and[5]for recent reviews). Unfortunately, 
lattice methods cannot be used in the real time formulation and because of this the 
study of quarkonium properties in hot QCD medium is a challenging problem. Still, 
our current understanding of QCD thermodynamics is largely based on lattice re- 
sults. In particular, our knowledge of the QCD equation of state, the nature of 
the transition and the corresponding transition temperature, and the response of 
the QCD medium to external static chromo-electric fields, all come from lattice 
QCD calculations. Most lattice calculations have been performed using so-called 
improved staggered fermion formulation. Variants of these formulations go under 
the acronyms asqtad, HISQ, p4, and stout (see Ref. [4] for a detailed discussion). 

As already mentioned in the introduction, hot strongly interacting matter un- 
dergoes a phase transition to a QGP, in which the dominant degrees of freedom 
are partonic. A rapid rise in the thermodynamic quantities (pressure and energy 
density), indicative of liberation of many new degrees of freedom, occurs around 
temperatures of 160 MeV and can be associated with deconfinement. As mentioned 
in the introduction, at temperatures of about 300 MeV the entropy density is only 
10 — 20% smaller than the entropy density of non-interacting gas of quarks and 
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gluons and these deviations from the ideal gas can be accounted for by using a re- 
summed weak coupling expansionPtHl The transition in QCD for physical light (u, 
d) and strange quark masses is found to be a crossover.^ Therefore, the definition 
of a transition temperature is non-trivial. As we mentioned in the introduction, 
the QGP is characterized by chiral symmetry restoration (broken in the vacuum 
and at low T). The chiral aspect of the QCD transition allows for a definition 
of a transition temperature since chiral symmetry is a good symmetry. The light 
quark masses are small compared to typical QCD scales, and therefore universality 
arguments allow one to relate the crossover temperature to the phase transition 
temperature in the limit of zero light quark masses P^Kfil It turns out that the phys- 
ical light quark masses are sufficiently small to allow such a definition and one 
obtains T c — 154(9) MeV This value is close to the inflection point of the renor- 
malized chiral condensate.^ However, chiral symmetry in the strange quark sector 
is restored at significantly higher temperatures. In addition to the equation of state, 
fluctuations of conserved charges have also been calculated in lattice QCD. At low 
temperatures, massive hadrons are the carriers of different conserved charges and, 
therefore, the fluctuations of conserved charges are suppressed by the corresponding 
Boltzmann factors. In fact, at low temperatures, fluctuations of conserved charges 
are well-described by a gas of hadrons and hadronic resonances EHE21 Similarly to 
the energy density and the entropy density, fluctuations of conserved charges show 
a rapid rise in a narrow temperature interval, and for T > 300 MeV are also well 
described by weak coupling techniques EUE1ED 

In the pure glue SU (N) theory there is a symmetry associated with the decon- 
finement transition, called the center of Z(N) symmetry. The partition function of 
this theory is invariant under gauge transformations, which leave the fields periodic 
up to a phase factor that is an clement of Z(N) group (this is a subgroup of the 
SU(N) gauge group called the center). This symmetry is spontaneously broken at 
high temperatures. Therefore, deconfinement is a true phase transition in SU(N) 
gauge theories. The corresponding phase transition temperature is 292(2) MeV 
(obtained through units of the string tension a is T/y/a = 0.629(3) assuming 
\[g = 465 MeV). The center symmetry is broken by the presence of dynamical 
quarks, and as we discuss in the next section, for physical values of the light quark 
masses it is not a good symmetry. Therefore, for QCD a deconfinement transition 
temperature, that is similar to the deconfining transition in pure gauge theory 
cannot be defined. 

2.2. Color Screening at High Temperatures 

The QGP is characterized by color screening. At high temperatures, color electric 
screening can be understood in terms of the low-momentum behavior of the tempo- 
ral component of the gluon propagator _Doo(fco,k). The corresponding self-energy, 
or the electric polarization function, in this limit is IIoo(0,k — > 0) = m,D with 
m,D = N/3gT-Jl + Nj/(2N) , mo being the leading order Debye massPHere Nf 
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is the number of light quark flavors: Nf = for pure gauge theory and Nf = 3 
for 2+1 flavor QCD (the effects of the strange quark mass are typically very small 
at high temperatures). The above result is gauge independent and, up to color and 
quark flavor factors, it is identical to the QED result, mry = eT/3P The trans- 
verse part of the static gluon propagator is not screened in perturbation theory, i.e. 
ilr(0,k — > 0) = . This, however, is not the case non-perturbatively, as we discuss 
this later. 

To study color electric screening non-perturbatively one uses Polyakov loops. 
As discussed in the previous subsection, the deconfming transition in SU(N) gauge 
theories is a true phase transition related to Z{N) symmetry. The order parameters 
of this phase transition are the expectation value of the Polyakov loop and the 
Polyakov loop correlator 

N T -1 

L(T) = (-W(f)), W(x) = [] U Q (r,x), (1) 

r=0 

C PL (r,T) = ±(TrW(r)TrW(0)). (2) 

The Polyakov loop transforms non-trivially under the Z(N) transformation. The 
expectation value of the Polyakov loop is zero in the confined phase and non-zero 
above the transition temperature. The qualitative change in the behavior of the 
Polyakov loop and its correlator above the phase transition temperature is related 
to color screening. The correlation function of the Polyakov loop is related to the free 
energy of a static quark anti-quark pair.^ But, as pointed out already in Ref. |2~T| 
the Polyakov loop and the Polyakov loop correlator require renormalization to be 
interpreted as the free energy of an isolated static quark or the free energy of a 
static quark anti-quark (QQ) pair. More precisely, these quantities are related to 
the free energy difference of a system with static QQ pair at some temperature and 
the system without QQ pair at the same temperature. Since in the zero temperature 
limit the free energy of a static quark anti-quark pair should coincide with the static 
potential, the Polyakov loop renormalization is determined by the normalization 
constant of the static potential, namely 

L ren {T) = exp(-c/(2T))L(T) = cxp(-F 00 (T)/(2T)) (3) 
C PL {r,T) =ex V (-F(r,T)/T + c/T), i^(T) = lim F(r,T), (4) 

r— >-oo 

where c is an additive normalization that ensures that the static potential has a 
certain value at a given distance. In the above expressions we made it explicit that 
the free energy of an isolated static quark, Fq , is one-half of the free energy of 
an infinitely-separated QQ pair. In the confined phase, the free energy of a static 
quark anti-quark pair is proportional to a(T)r at large distances r, as expected in 
QCD. The effective temperature-dependent string tension <r(T) is non-zero below 
the phase transition temperature!^^ Consequently, the free energy of an isolated 
static quark is infinite and L ren (T) = 0. In the deconfined phase, the free energy 
of an infinitely separated QQ pair, F^IT) , is finite due to color screening. In 
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particular, at high temperatures Foo = — (N 2 — l)a s mD / '(2JV) at leading order. 
This result follows from the infrared behavior of the self energy n 00 (0, k) discussed 
above. The next-to-leading order correction for has been calculated in Refs. l24l 
and [251 and was found to be small. For the free energy of a static QQ pair at distance 
r we have at leading order^ 

la 2 N 2 -l 
F(r, T ) = ~-fp-^ exp(-2m D r) ^—a s m D . (5) 

Note, that this result is in contrast to the free energy of static charges in QED 
which is — aexp(— m£>r)/r — amo, and is the consequence of non-Abelian nature 
of interactions. A static quark anti-quark pair could be either in a singlet or in an 
octet state, and thus (in a fixed gauge) one can define the so-called color singlet 
and octet free energjJ2M22l@ 

exp(-*i(r,T)/T + c/T) = 1- (Tr[W^x)W (y)}) , (6) 

exp(-F 8 (r,T)/r + C /T) = -J—^w^x) TrW(y)) 

- N{ ^_ 1) (^[W\x)W(y)]). (7) 

One can view the above correlators as correlation functions of singlet and octet 
static meson operators of size r in Coulomb gauge evaluated at r = l/TpSE^Then, 
the Polyakov loop correlator can be written as the thermal average over the singlet 
and the octet contributions JUIHsMH 

1 N 2 - 1 

CpL(r,T) = _ex P (-F 1 (r,T)/T + c/T) + -^ 5 - ffl q)(-F 8 (r,T)/r + c/T). (8) 

The singlet and octet free energies can be calculated at high temperature in leading 
order hard-thermal-loop (HTL) approximation^ resulting in 

jpf t\ N 2 -l a s (N 2 - l)a s m D 

Fi{r,T) = ———exp(-m D r) — , (9) 

Fs(r,T) - — — cxp(-m D r) — . (10) 

The singlet and octet free energies are gauge independent at this order. Inserting 
the above expressions into Eq. and expanding in a s , one can see that the 
leading order contributions from the singlet and octet channels cancel each other 
and one recovers the leading order result for F(r,T) in Eq. ([5]). Thus the non- 
Coulombic nature of the QQ free energy is due to the partial cancellation of the 
singlet and octet contributions at leading order. For a long time it was not clear 



a The terms color singlet free energy and octet free energy are misleading as only F(r, T) has 
the meaning of the free energy of static quark anti-quark pair. Furthermore, at low temperatures 
both Fi and Fg are determined by color singlet asymptotic statesPH The color screening of the 
static quark quark interaction was discussed in Ref. 1291 while color screening of charges in higher 
representations was studied in Refs. [30l and l3ll 
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Fig. 2. The Polyakov loop as function of temperature in 2+1 flavor QCrJS*E^and in pure gauge 
theory The right panel shows the free energy of a static quark Fq = Fac/2 as function of 
temperature and compared with the hadron resonance gas model of static-light hadrons^l 



how to generalize the decomposition of the free energy into the singlet and octet 
contributions beyond leading order. Recently, using the effective theory approach 
that will be discussed below, it was shown that this decomposition indeed holds 
at short distances.^ As we will see below, the singlet free energy turns out to be 
rather useful when studying screening numerically in lattice calculations. 

As already mentioned, the QCD partition function does not possesses Z(3) sym- 
metry. This symmetry is broken by the quark contribution. In physical terms this 
means that a static quark anti-quark pair can be screened already in the vacuum 
by light dynamical quarks. The free energy of a static QQ pair does not rise lin- 
early with distance r indefinitely but saturates at some distance, i.e. we see string 
breaking. For light dynamical quarks the corresponding free energy, i 7 ^ , at very 
low temperatures is twice the binding energy of static-light meson and thus is pro- 
portional to Aqcd- As the temperature increases excited static light mesons, as 
well as baryon states with one static quark will contribute to the expectation value 
of the Polyakov loopP^ Using the spectrum of known hadron states with a heavy 
quark or using the spectrum of static light hadrons, the expectation value of the 
Polyakov loop can be calculated in the approximation of a non-interacting gas of 
these hadronsP 1 ^ In Fig. [2] we show the recent lattice results for the renormalized 
Polyakov loop in 2+1 flavor QCD with physical quark masses extrapolated to the 
continuum IimhPEl and compare these to the results in pure gluodynamics.^^ 
The temperature scale in pure gauge theory was set using the value yfa = 465 MeV 
for the zero temperature string tension. At high temperatures the behavior of the 
Polyakov loop in pure glue theory and QCD look similar, at least qualitatively. 
The difference appears when approaching the transition region. The Polyakov loop 
behaves smoothly across the transition region. We also show the comparison of the 
free energy of a static quark, Fq — Foe/2 , with the prediction of a resonance gas of 
static-light states,^ which seems to work well within the uncertainties of the model 
up to temperatures of about 140 MeV. 
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The free energy of a QQ pair, as well as the singlet free energy have been studied 
in 2+1 flavor QCD using the HISQ action and 16 3 x 4 and 24 3 x 6 lattices.^ The 
numerical results for the singlet free energy are shown in Fig. [3] At short distances 
the singlet free energy agrees with the zero temperature potential calculated in 
Ref. [TBI while at large distances it approaches a constant value Foo (T) , equal to 
the excess free energy of two isolated static quarks. As the temperature increases, 
the deviation from the zero-tcmpcrature potential shows up at shorter and shorter 
distances as the consequence of color screening. For temperatures T > 200 MeV 
there are large deviations from the zero temperature potential already at distances 
r that are similar to charmonium size. Thus we would expect strong modification 
and/or dissolution of charmonium states. To explore the screening behavior in Fig. 
[3] we also show the combination S(r, T) = r ■ (F(r, T) — F^T)) which we call the 
screening function. The screening function should decay exponentially. We indeed 
observe the exponential decay of this quantity at distances larger than 0.8/T . Thus 
at high temperatures the behavior of the singlet free energy expected from the weak- 
coupling calculations seems to be confirmed by lattice QCD, at least qualitatively. 
Let us also mention that at high temperatures the behavior of the singlet free energy 
is similar to that observed in pure gauge theorypjEHD 

In Fig. 0] we show lattice results for the free energy of a static QQ pair as 
function of the distance at different temperatures. At short distances and low tem- 
peratures the free energy is expected to be dominated by the singlet contribution 
and we expect it to be equal to the zero temperature potential up to the term 
T In 9 , term coming from the normalization (see the discussion in the previous sec- 
tion). Therefore, in the figure, the numerical results have been shifted by —Tin 9 . 
Indeed, for the smallest temperature and the shortest distances F(r,T) — Tin 9 
is equal to the zero temperature potential shown as the dashed black line. At 
higher temperatures F(r,T) is very different from the zero-temperature poten- 
tial. At large distances the free energy approaches a constant value Foo(T) that 
decreases with increasing temperatures, confirming expectations (see discussions 
above). The temperature-dependence of F(r,T) is much larger than that of the 
singlet free energy. This is presumably due to the partial cancellation of the sin- 
glet and octet contribution discussed above. To verify this assertion we calculated 
F(r,T) — -F'oo(T) using the numerical data for F\{r,T) — F^T) and the leading 
order relation (F\ (r, T) - F 00 (T))/(F 8 (r, T) - F^ (T)) = -8 . The corresponding 
results are shown in the right panel of Fig. |4] As one can see from the figure, the 
numerical data for F(r, T) are in reasonable agreement with the ones reconstructed 
from this procedure. The reconstruction works better with increasing temperature. 
Thus the expected cancellation of the singlet and octet contributions to the free 
energy of static QQ pair seems to be confirmed by lattice calculations. The numer- 
ical results discussed above are similar to the one obtained with p4 actiorP^"^ for 
T > 200 MeV. 

The qualitative features of the singlet free energy discussed above are not specific 
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Fig. 4. The free energy of static QQ pair (left) and the difference F(r, T) — Foo{T) (right) calcu- 
lated with HISQ action as function of quark separation distance r at different temperatures. In the 
right panel the filled symbols correspond to the lattice data, while the open symbols correspond 
to the values reconstructed from the singlet free energyPSl The legend in the left panel is the same 
as in Fig. [3] (left). 

to the Coulomb gauge. In fact, the singlet free energy defined in terms of Wilson 
loops shows very similar behavior to the one calculated in the Coulomb gauged 
Furthermore, the Debye screening mass can also be defined from the long distance 
behavior of the electric gluon propagator PHHII] Calculations in 517(2) and 517(3) 
gauge theories show no gauge dependence of the extracted Debye mass within sta- 
tistical errors P3EU The extracted screening masses are in agreement with the ones 
obtained from the singlet free energy. Color magnetic screening can be also studied 
in this approachp^H^ though the long distance behavior of the static transverse 
gluon propagator is not always exponential^"^ and thus the corresponding screen- 
ing mass cannot be defined. 



2.3. Quarkonium Spectral Functions and Euclidean Correlators on 
the Lattice 

In-medium meson properties are encoded in meson spectral functions. The spectral 
function a(cu,p) for a given meson channel H in a system at temperature T can 
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Table 1. Meson states in different channels for light, charm, and bottom quarks. 
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be defined through the Fourier-transform of the real time two-point function D > 
and D < of the meson current Jh , or equivalently, as the imaginary part of the 
Fourier-transformed retarded correlation functional 

a(cj,f) = ^-(D>(cj,f) - D<(u>,0)) 



2tt 
1 

7T 

oc 



ImD H (cu,p) 
dt 



d 3 xe ^t-tp-x D >(<) 



{t,x) 



{J H {t,x)J H (0,0)) 
{J H (0,(})J H (t,x)),t >0 



D>(t,x) 
!><(«, i?) 

Here we study local meson operators of the form 

J H {t,x) = q(t,x)T H q(t,x) 
where q(t, x) is the quark field operator and 



(11) 
(12) 
(13) 
(14) 



for scalar, pseudo-scalar, vector, axial-vector and tensor channels. The relation of 
these quantum number channels to different meson states is given in Table [1] 

The correlators D^ < \t,x) satisfy the well-known Kubo-Martin-Schwinger 
(KMS) condition^ 

D>{t,x)=D<{t + i/T,x). (15) 
Inserting a complete set of states and using Eq. (|15[) . one obtains the expansion 

^P) = iE m , n e- E " /T * 
\(n\ J H (0)|m}| 2 (<5 4 (^ + fc» - fc™) - ^ + fc™ - *£)) (16) 

where Z is the partition function, /c"( m ) refers to the four-momenta of the state 
\n(m)) and p M = (u),p). 

A stable meson state contributes a 5-function-like peak to the spectral function: 

cr(u,p) = \(Q\J H \H)\ 2 e(cj)6(p 2 - M 2 ), (17) 

where Mh is the mass of the state and e(po) is the sign function. For a quasi-particle 
in the medium, one obtains a smeared peak, with the width being the thermal 
width. As one increases the temperature, the width increases and at sufficiently 
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high temperatures, the contribution from the meson state in the spectral function 
may be sufficiently broad, and thus it is no longer meaningful to speak of it as a 
well-defined state. 

In finite temperature lattice calculations, one calculates Euclidean time propa- 
gators, usually projected to a given spatial momentum: 



G(r,p)= / d 6 xe* p - x {J H (T,x)J H (0,0)}. (18) 

This quantity is an analytical continuation of D^(t,p) 

G(r,p) = D>(-ir,p). (19) 

Using this together with the periodicity of the time direction, we obtain the following 
integral representation for the Euclidean time correlator 



G{t,p) = I dh)a{uj,p)K(u),T), 
Jo 



coshMr - 1/2T)) 
K ^ T)= sinter) ■ (20) 

This equation is the basic equation for extracting the spectral function from meson 
correlators. Equation (|20[) is valid in the continuum. Formally, the same spectral 
representation can be written for the Euclidean correlator calculated on the lat- 
tice G lat (T,p). The corresponding spectral function, however, will be distorted by 
the effect of the finite lattice spacing, in particular, the spectral function is zero 
above certain energy uj > Lo max . These distortions have been calculated in the free 
theory™] 

Early attempts to reconstruct charmonium spectral functions from Euclidean 
quarkonium correlators in quenched approximation (pure glue theory) using the 
Maximum Entropy Method (MEM) have been presented in Refs.[49l[50j and[5TJ The 
reconstructed spectral functions contained some peak-like structures that survived 
up to temperatures as high as 1.6 T c and even higher, to 2T C . These structures have 
been interpreted as the charmonium ground state surviving into the deconfined 
phase. This is contrary to the picture of color screening obtained from the singlet 
free energy. Similar results have been presented in Ref. [52] for two flavor QCD. 
It turned out, however, that the above conclusion was premature: As we discuss 
below, quarkonium correlators are not particularly sensitive to modifications in the 
spectral function due to the dissolution of the quarkonium bound states. Detailed 
systematic studies revealed that with more realistic priors in the MEM analysis 
there is no indication of the existence of the bound states^"^ (for the discussion of 
MEM systematics also see Ref. 55). 



2.3.1. Temperature- dependence of the Euclidean Correlators 

Since the analysis of the quarkonium spectral functions using MEM turns out to 
be quite complicated, it is important to understand the temperature-dependence 
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of the quarkonium correlation functions, as well as try to identify possible sources 
of the temperature-dependence that are related to melting of the bound states. 
From Eq. (|20[) it is clear that the temperature-dependence of the meson correla- 
tion functions comes from two sources: the trivial temperature-dependence of the 
integration kernel K (r, uS) and the temperature-dependence of the spectral function 
(t(uj, T). To get rid of the first trivial temperature-dependence, one can consider the 
reconstructed correlation function 

/>oo 

G rec (r,T)= / (L)a(uj,T=0)K(uj,T). (21) 
Jo 

If the spectral function does not change across the deconfinement transition 
G(T,T)/G rec (T,T) should be unity. Deviations of this ratio from unity indicate 
a temperature-dependence of the spectral functions. The ratio G(r, T)/G rec (r, T) 
was first studied in Ref. 02] and subsequently in Refs. [551 [S7J [551 an d E31 It was 
found that in the pseudoscalar channel this ratio stays close to one, showing only 
small temperature-dependence. This result seemed to support the early conclusion 
based on the spectral functions extracted from MEM, that ground state charmo- 
nium survives in the deconfined phase up to temperatures 1.6 T c . In the scalar and 
axial- vector channels large temperature dependence in G(r 1 T)/G rec (r 1 T) was ob- 
served. This was interpreted as melting of IP charmonium states and fit well into the 
expected picture of sequential melting. However, as already mentioned above, this 
conclusion was premature. In the deconfined phase there is an additional contribu- 
tion to the spectral functions at very low frequency. This can be seen by calculating 
the spectral function in the free theory, where a contribution proportional to lj5{uj) 
appears in all but the pseudoscalar channel.^ This is addition to the continuum 
contributions which start at u) = 2m q , where m q is the heavy quark mass. The 
delta function gets smeared once the interactions of the heavy quarks with the 
medium are taken into account. In particular, in the vector channel, where this low 
u) structure is related to the heavy quark diffusion we have^ 

uS(u) --5^-2 , r, = T/(Dm q ), (22) 

D being the heavy quark diffusion constant. The schematic structure of the spec- 
tral functions is shown in Fig. [5j The area under the transport peak at low uj is 
given by the quark number susceptibility \q- F° r large quark mass the transport 
peak is narrow and is well-separated from the high energy part of the spectral func- 
tion, which corresponds to bound states and/or unbound heavy quark anti-quark 
pairs (continuum) . One expects a similar structure in the spectral function in other 
channels as well. Thus we can write 

cr(w,T) =cri ow {u),T) + a high (w,T), (23) 
G(r, T) = G low (r, T) + G high (r, T) . (24) 

Since the peak at low u is narrow, we expect that the derivative of Gi ow (r, T) 
with respect of r is small, G' 1ow (t,T) ~ 0. Lattice calculations in quenched QCD 
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Fig. 5. The schematic structure of the vector spectral function for heavy quarks. The arrow 
indicates the value of the spectral function at zero frequency. 



show that the ratios G 1 (r) / G' rec (T) have no strong temperature-dependence and are 
close to one for all temperatures, including the highest temperature of 3T c P^This 
means that the large temperature-dependence seen in G/G rec for the scalar and 
axial- vector channels comes mostly from the temperature-dependence of Gi ow (r, T) 
and is not related to melting of P-wave charmonium states. At the same time 
Ghigh, which is sensitive to the bound state properties, does not show significant 
temperature-dependence for all of the channels. A question thus arises, as to how 
the observed temperature- (in)dependence of Gh%gh(j,T) is related to the expected 
in- medium modification of the heavy quark potential and the melting of quarkonium 
states at sufficiently high temperatures. This question is addressed below in context 
of potential models in Sections 12.4.21 and 12.4.31 



2.3.2. Spatial Meson Correlators 

In lattice QCD one also calculates a meson correlation function in one of the spatial 
directions, say z 

G(z,T) = j dxdy jj dr{J(x, y, z, r) J(0, 0, 0, 0)). (25) 

The spatial correlation function is related to the meson spectral function at non-zero 
spatial momentum 

G(z,T)= [°° d Pz e^ z [°° du aiuJlPz ' T) . (26) 



Thus the temperature-dependence of the spatial correlation function also provides 
information about the temperature-dependence of the spectral function. Unlike the 
temporal correlators, spatial correlators can be calculated for arbitrarily large sep- 
aration and, therefore, can be more sensitive to in-medium modification of quarko- 
nium spectral functions. Medium effects are expected to be the largest at distances 
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which are larger than 1/T. At these distances G(z, T) decays exponentially. This ex- 
ponential decay is governed by the so-called screening mass, M scr . If there is a lowest 
lying meson state of mass M, i.e. the spectral function can be well-approximated 
by Eq. (fTTj) . then the long-distance behavior of the spatial meson correlation func- 
tion is determined by the meson mass, i.e. M scr — M . At very high temperatures, 
the quark and anti-quark are not bound and the meson screening mass is given 
by 2^J(irT) 2 + m 2 , where m q is the quark mass and irT is the lowest Matsubara 
frequency. Therefore, a detailed study of spatial meson correlators and screening 
masses can provide some information about the melting of meson states at high tem- 
peratures. Spatial meson correlators and screening masses have been calculated in 
the light hadron sectoi^ and for charmoniumP*2l The changes in the correlators and 
screening masses around the transition temperature are much smaller for quarkonia 
than for the light mesons. This fits into the picture that ground state of charmonium 
does not melt at the crossover temperature. However, for T > 200 MeV we see large 
change in spatial charmonium correlators at distances zT > 1/2 . Furthermore, for 
T > 300 MeV the behavior of the screening masses and their dependence on the 
spatial boundary conditions is compatible with the picture of unbound quarksP*2l 
This corroborates the discussion in the previous subsection that there is no evidence 
for survival of charmonium states at temperatures of 1.6 T c . 



2.4. Effective Field Theory Approach 



There are different scales in the heavy quark bound state problem related to the 
heavy quark mass m q , the inverse size ~ m q v ~ 1/r , and the binding energy 
m q v 2 ~ a s /r . Here v is the typical heavy quark velocity in the bound state and 
is considered to be a small parameter. Therefore, it is possible to derive a sequence 
of effective field theories using this separation of scales (see Refs. [63] and [64] for 
recent reviews). Integrating out modes at the highest energy scale ~ m q leads to an 
effective field theory called non-relativistic QCD (NRQCD), where the pair creation 
of heavy quarks is suppressed by powers of the inverse mass and the heavy quarks 
are described by non-relativistic Pauli spinorsP^ At the next step, when the large 
scale related to the inverse size is integrated out, the potential NRQCD or pNRQCD 
appears. In this effective theory the dynamical fields include the singlet S(r,R) and 
octet 0(r, R) fields corresponding to the heavy quark anti-quark pair in singlet and 
octet states respectively, as well as light quarks and gluon fields at the lowest scale 



The Lagrangian of this effective field theory has the form 
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4 ^ 



^2 Ql iJpi: 
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(27) 
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Here the dots correspond to terms which are higher order in the multipole expan- 
sion.^ The relative distance r between the heavy quark and anti-quark plays a 
role of a label, the light quark and gluon fields depend only on the center-of-mass 
coordinate R. The singlet V s (r) and octet V a (r) heavy quark potentials appear as 
matching coefficients in the Lagrangian of the effective field theory and, therefore, 
can be rigorously defined in QCD at any order of the perturbative expansion. At 
leading order 

V s (r) = -^, V (r) = l^ (28) 
6 r or 

and Va — Vb — 1- One can generalize this approach to finite temperature. How- 
ever, the presence of additional scales makes the analysis more complicated.^ The 
effective Lagrangian will have the same form as above, but the matching coeffi- 
cients may be temperature-dependent. In the weak coupling regime there are three 
different thermal scales: T, gT and g 2 T. The calculations of the matching coeffi- 
cients depend on the relation of these thermal scales to the heavy quark bound-state 
scales!^ To simplify the analysis the static approximation has been used, in which 
case the scale m q v is replaced by the inverse distance 1/r between the static quark 
and anti-quark. The binding energy in the static limit becomes V — V s ~ Na s /(2r). 
When the binding energy is larger than the temperature the derivation of pNRQCD 
proceeds in the same way as at zero temperature and there are no medium mod- 
ifications of the heavy quark potentialPS But bound state properties will be af- 
fected by the medium through interactions with ultra-soft gluons. In particular, 
the binding energy will be reduced and a finite thermal width will appear due to 
medium-induced singlet-octet transitions arising from the dipole interactions in the 
pNRQCD LagrangiarP*] (c.f. Eq. (|27p). When the binding energy is smaller than one 
of the thermal scales, the singlet and octet potential will be temperature-dependent 
and will acquire an imaginary partP^ The imaginary part of the potential arises 
because of the singlet-octet transitions induced by the dipole vertex as well as due 
to the Landau damping in the plasma, i.e. scattering of the gluons with space-like 
momentum off the thermal excitations in the plasma. In general, the thermal cor- 
rections to the potential go like (rT) n and (m£>r) n p3 where m_o denotes the Debye 
mass. Only for distances r > 1/mc there is an exponential screening. In this region 
the singlet potential has a simple form 

V(r)=--^-e- mDr + i-a T - [°° dx sin ( m ^ rx ) _ 1 a ( mn+lT ) 

Vs[) ~ 3r 6 + 3 s rm D J X (x 2 + l) 2 3 <*'V" D +tl), 

The real part of the singlet potential coincides with the leading order result for 
the free energy The imaginary part of the singlet potential in this limit was first 
calculated in Ref. 1671 For small distances, the imaginary part vanishes, while at large 
distances it is twice the damping rate of a heavy quark!^! This fact was first noted 
in Ref. |69] for thermal QED. The effective field theory at finite temperature has 
been derived in the weak-coupling regime assuming the separation of the different 
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thermal scales and also Aqqd- In practice, the separation of these scales is not 
evident and one needs lattice techniques to test this approach. 

2.4.1. Weak Coupling Results on Quarkonium Properties 

The binding energies and the thermal width of quarkonium have been calculated 
using pNRQCD in the weak coupling regime, assuming the following hierarchy of 
scales: m q a s »T» m q a 2 s 3> mj)P This regime may be of relevance for the T, 
the ground state bottomonium in the QGP. It was found that the thermal shift 
in the binding energy is proportional to T 2 /(a s m g ), while the thermal width is 
proportional to a s T. In Ref. 71 it has been shown that the thermal width obtained 
in effective field theory framework naturally incorporates previous estimates of the 
width that rely on kinetic theory and the gluo-dissociation cross-section.^ The 
structure of the spin-dependent correction to the finite temperature potential has 
been studied in the weak coupling regime,^ as well as the thermal correction to the 
potential for bound states that move with respect to the rest frame of the plasmaP^ 
In particular, it was found that above certain temperature the Landau-damping is 
not effective as a mechanism for bound state dissociation.^ 



2.4.2. Beyond Weak Coupling: Potential Models 

At increasingly high temperates quarkonium binding energies are significantly re- 
duced from their vacuum values and eventually vanish. This is called the point of 
zero binding. In this situation the binding energy is the smallest scale in the prob- 
lem and all the medium effects can be incorporated in the potential. If we start 
from pNRQCD and neglect the octet-singlet interactions the dynamics of the sin- 
glet quark anti-quark field is given by the free field equation, which is the form of 
the Schrodinger equation (e.g. see discussions in Ref. [S3]). Therefore, the calculation 
of the correlation function of singlet fields in real time, or equivalently the calcu- 
lation of quarkonium spectral function in the non-relativistic limit, is reduced to 
calculating the Greens-function of the Schrodinger equation^ 



V 2 + V(r,T) - E 

m„ 



G nr (f,P,E,T) =S 3 (r-r') 



a(u,T) = -lmG nr (r,r',E)\ p _ r -,_ Q E = oj-2m. (29) 

The potential entering the above equation is complexp^^ and, as discussed above, 
its form is quite complicated even in the weakly coupled regime, since the form of 
the thermal corrections depends on the relation of the scales 1/r, T, mj ~ gT. 
Furthermore, in the relevant temperature regime the separation of the above scales 
does not hold and the potential is effected by the non-perturbative scales, such as 
g 2 T and Aqcd ■ Therefore, in this case we have to rely on the lattice calculations 
to constrain the form of the potential. 
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Fig. 6. The S-wave charmonium (left) and bottomonium (right) spectral functions calculated 
for quenched QCD using a lattice QCD inspired potential: The insets show the ratio G/G rec 
obtained from the calculated spectral functions and compared to lattice QCD data. 



2.4.3. Potential Models and Euclidean Correlation Functions 

A simple potential model based on lattice QCD data on the singlet free energy has 
been proposed to study quarkonium spectral functions and explain the observed 
temperature-independence of the quarkonium correlators! 76 ^ 77 ! The imaginary part 
of the potential has been neglected in this work, while the real part of the potential 
was based on lattice QCD data on the singlet free energy. In this model all quarko- 
nium states except the ground state bottomonium melt in the QGP. However, the 
interaction between the heavy quark and anti-quark leads to large enhancement 
of the spectral function in the threshold region. This is shown in Fig. [6] for the 
case of S-wave charmonium spectral function in quenched QCD. One can see, that 
there is no bound state already at 1.2 T c , while the correlation functions do not 
show large temperature-dependence. The ratio G/G rec is close to unity, in agree- 
ment with the lattice findings. Similar results have been found also for the P-wave 
quarkonium EsES This result is understood in that the large threshold enhance- 
ment compensates for the absence of bound states and therefore the correlators do 
not change much. Even though the difference between the vacuum and in-medium 
spectral function increases with increasing temperatures, the Euclidean time extent 
1/(2T) becomes smaller. Therefore, it is difficult to see this change in Euclidean 
time correlators. Similar results have been reported also in Ref. [78] Furthermore, 
the above findings do not depend on the precise choice of the potential, provided it 
is compatible with lattice QCD findings at short distanced 7 - 6 -" 7 - 7 -' (see also discussions 
in Ref. I79|) . This analysis resolved the apparent puzzle between the strong modi- 
fication of the potential and very small temperature dependence of the Euclidean 
quarkonium correlators. 
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2.4.4. Static Energy in Lattice QCD 

As already mentioned, in the limit of small binding energy all the thermal scales, 
as well as Aqqd , can be integrated out, and the potential can be approximated by 
the static energy. The static energy has been calculated in perturbation theory^ 
and is complex-valued. Since the validity of perturbation theory is not evident for 
the relevant temperature range, it is interesting to calculate the static energy on 
the lattice. In the previous section we considered the correlation function of a static 
QQ pair evaluated at Euclidean time t = 1/T. These correlators are related to the 
free energy of a static QQ pair. One can consider Wilson loops for Euclidean times 
t < 1/T, which have no obvious relation to the free energy of a static QQ pair. 
Wilson loops at non-zero temperature have been first studied in Refs. [80] and [81] 
in connection with the heavy quark potential at non-zero temperature. A spectral 
decomposition for the Wilson loops has been conjectured: 

/>oo 

W(t,t)= a(u,r,T)e- UT . (30) 
Jo 

At zero temperature, the spectral function is proportional to a sum of delta func- 
tions a(r,uj) — ^2 n c n S(E n (r) — w) . At high temperatures, the spectral function is 
proportional to a sum of smeared delta functions. The position and width of the 
lowest peak in the spectral function are related to real and imaginary part of the 
potential, respectivelyP^ED Motivated by this, Wilson loops have been calculated on 
finite temperature lattices in 2+1 flavor QCD using the HISQ action with physical 
strange quark mass and light quark masses that correspond to a pion mass of 160 
MeV in the continuum limit ISZllSS] p rom the analysis of the Wilson loops, the real 
part of the static energy of QQ pair has been extracted. Unfortunately, the Wilson 
loops are not very sensitive to the imaginary part of the static energy. The corre- 
sponding results at two temperatures are shown in Fig. [71 The real part of the static 
energy is larger than the singlet free energy. As one can see from the figure, screen- 
ing effects are not present at T = 178 MeV but can be clearly seen for 220 MeV. 
This finding is in accord with the expectations from the temperature-dependence 
of Foo(T) which seems to be described well by hadronic degrees of freedom for 
T < 200 MeV (see the previous section). 

2.4.5. Spectral Functions in a Realistic Potential Model 

As discussed above, the in-medium corrections to the potential have both real and 
imaginary parts. For very small couplings the imaginary part of the potential is 
actually the dominant source of bound state dissolution. Therefore, any realistic 
calculation of quarkonium spectral functions should include both real and imaginary 
parts for the temperature-dependent potential. The hierarchy of scales assumed in 
weak coupling calculations may not be satisfied and non-perturbative effects may 
play an important role (see previous section). Thus, one has to rely on the lattice 
results. But, as mentioned before, the imaginary part of the potential cannot be 
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extracted from present lattice calculations. Therefore, in Ref.[83lthe imaginary part 
of the potential was evaluated according to Eq. (p9|) using the prescription for the 
temperature-dependent running coupling constant from Rcf. [731 For the real part of 
the potential lattice information on color screening was used. Namely, the potential 
was exponentially screened for distances r > r scr — 0.8/T, and at smaller distances 
this was taken to be equal to the zero temperature potential. This construction 
gives the maximal real part of the potential that still compatible with the data. 
The real and imaginary part of the potential used in this analysis is shown in Fig. 
[H We also show the static energy calculated on the lattice for the matching values 
of the temperature. Interestingly, the real part of the potential is close to the lattice 
data on the static energy. The calculated S-wave charmonium and bottomonium 
spectral functions are shown in Fig. [9l As one can see from the figure, all of the 
quarkonium states, except the IS bottomonium, are dissolved for T > 250 MeV. 
The T(liS') may survive up to temperatures of 450 MeV and its melting will only 
occur at higher temperatures. We would like to re-iterate, that the key ingredient in 
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Fig. 9. Charmonium (left) and bottomonium (right) spectral functions calculated in a potential 
model with complex potentialJSH 



the observed melting pattern is the presence of the imaginary part of the potential; 
neglecting it would lead to larger dissociation temperatures (see e.g. Ref. |84|). The 
above estimates of the melting temperatures are in agreement with earlier estimates 
presented in Ref. [77] (note, that T c = 204 MeV in those pure gauge calculations). 
Also, the above spectral functions agree qualitatively with the ones found in a T- 
matrix approach that also includes some of the effects of the imaginary potential.^ 



2.4.6. Lattice NRQCD at Non-zero Temperature 

For the most interesting cases the heavy quark mass is well-separated from other 
scales. Therefore, NRQCD is always a good effective theory and can be combined 
with non-perturbative approaches like lattice QCD. A lattice version of NRQCD has 
been successfully used to study bottomonium spectral functions non-perturbatively 
(see e.gP^). This approach can also be used to study bottomonium spectral func- 
tions at non-zero temperatures.^^ Studying bottomonium spectral functions using 
relativistic heavy quarks turns out to be difficult because of the noise, which is not a 
problem in the non-relativistic formulation. In the NRQCD approach heavy quarks 
are not subject to anti-periodic conditions and, therefore, the meson correlator can 
be studied up to twice larger distances r < 1/T and has a simpler spectral decom- 
position 

/oo 
duja{uj,T)e~ UT . (31) 
-oo 

Due to the fact that gauge fields are still periodic, the spectral function can be non- 
zero for negative energies.^ In this approach there is no zero-mode contribution. 
This, and the fact that correlators can be studied at larger distances, makes the 
correlators more sensitive to the medium-modification of the spectral functions. The 
temperature dependence of P-wave bottomonium correlators appears to be strong 
at temperatures of about 280 MeV, which is indicative of melting of P states.^ At 
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a temperature of 400 MeV, the behavior of the P-wave correlators is consistent with 
propagation of unbound quarks.^ In contrast, the S-wave bottomonium correlators 
do not show significant temperature dependence which is consistent with survival 
of bottomonium states to temperatures T > 400 MeV. 

3. Quarkonium Phenomenology 

3.1. The Anisotropic (Viscous) Quark Gluon Plasma 

State of the art dynamical models of QGP spatiotemporal evolution rely on some 
form of viscous hydrodynamical treatment. Over a decade ago it was shown that 
ideal relativistic hydrodynamics was able to reproduce the soft collective flow of the 
matter and single particle spectra produced in relativistic heavy ion collisions jSsHUl 
However, one expects on general grounds that there is a lower bound on the QGP 
shear viscosityP^J] Based on this realization, there was subsequently a concerted 
effort to develop a systematic framework for describing the soft dynamics of the 
QGP using relativistic viscous (non-ideal) hydrodynamics p03tl22l 

Traditional viscous hydrodynamics approaches assume that the system is close 
to thermal equilibrium and hence is also very close to being isotropic in momentum 
space. However, one finds during the application of these methods that this assump- 
tion breaks down at the earliest times after the initial impact of the two nuclei, due 
to large momentum-space anisotropies in the pt~Pl plane, which can persist for 
the entire lifetime of the plasma.^SiHinU The level of the observed momentum-space 
anisotropy increases as the assumed shear viscosity of the QGP increases. Similar 
conclusions have been obtained in the context of strongly coupled systems, where it 
has been shown that one achieves viscous hydrodynamical behavior at early times, 
however, the system still possesses large momentum-space anisotropies which per- 
sist throughout the evolution.EzKH] As a result, it is necessary to include the effect 
of momentum-space anisotropies on the screening of heavy quarkonium. 

This is accomplished minimally by introducing a momentum-space anisotropy 
parameter, £, which is related to the ellipticity of the distribution function in the 
Pt~Pl plane ElllliSI m t the underlying one-particle partonic distribution functions. 
The state of the QGP is then described by a local average momentum scale, A(t, x) 
(transverse temperature), and a local momentum anisotropy parameter £(t, x). If 
the system is close to being isotropic in momentum space, then £ can be mapped to 
the plasma shear, n, and 2nd order viscous hydrodynamics is recovered; however, 
if the system is highly anisotropic it is possible take moments of the Boltzmann 
equation and derive a set of coupled partial differential equations! 105 ^ 106 ^ 119 " 124 ! The 
resulting dynamical equations allow one to vary the assumed value of the QGP shear 
viscosity in order to assess its impact on the dynamical evolution of the system. 

Additionally, one must revisit the calculation of the heavy quark potential in 
the presence of momentum-space anisotropies. This was first considered in RefP^ 
in which the real part of the heavy quark potential was calculated. Shortly af- 
terwards, the imaginary part of the heavy quark potential in the the presence of 
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Fig. 10. Real and imaginary parts of the binding energy of the T(ls) (left) and T(2s) (right) 
states as a function of the temperature for a momentum-space isotropic QGP. Dashed vertical 
lines indicate the temperature at which the real and imaginary parts are equal. Numerical results 
taken from Ref. 11301 



momentum-space anisotropics was calculatecP^SHHIl anc j phenomenological poten- 
tial models incorporating this effect were constructed P^ESD The chief effect of 
momentum-space anisotropies is to reduce Debye screening in the plasma@ This in 
turn has the effect that heavy quarkonium states can survive to higher tempera- 
tures in the presence of momentum-space anisotropies. The effect of mild (£ = 1) 
momentum space anisotropies is indicated in the right column of Table [H From 
this we can see that, for all states shown, the dissociation temperature increases as 
the level of momentum-space anisotropy increases. We consider a state dissociated 
when its width is larger than its binding energy, because a state then decays faster 
than it binds. The temperatures shown in [5] are the temperatures where thermal 
medium effects of ~ T are comparable with the binding energy. This condition at 
finite temperatures is more realistic than the zero-binding condition, often quoted 
in the literature, which is sensible only at zero temperature. 



3.2. A Phenomenological Potential for an Anisotropic Plasma 

In Ref. 11301 a complex- valued potential model was constructed in order to address 
the question of T suppression, in particular. The imaginary contribution to the 
potential was taken to be given by the leading order perturbative result jS3ll25|i26| 
The real part was taken to be given by either the free energy or internal energy, with 
the free energy, F, given by a Karsch-Mehr-Satz potential forrrP^ and the internal 
energy given by U = F — TS. In both cases, lattice results were used to fix empirical 
constants such as the string tension. For simplicity, herein we focus on the results 
obtained using the more binding internal energy model, since the free-energy-based 
model seems to give too much suppression at both RHIC and LHC energies. 

b Another interesting, but subleading, effect is that momentum-space anisotropies cause a splitting 
between £ ^ states. 
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Table 2. Estimates of the isotropic and anisotropic dissociation scales for the 
J/1/1, Xcli f(l s ); T (2s), T(3s), Xbl> an d Xb2- Estimates are taken from Refs. lT29l 

[T501 



State 


Isotropic QGP (5=0) 


Anisotropic QGP (5=1) 




307 MeV 


374 MeV 


Xcl 


< 192 MeV 


210 MeV 


T(l«) 


593 MeV 


735 MeV 


T(2s) 


228 MeV 


290 MeV 


T(3 S ) 


< 192 MeV 


< 192 MeV 


Xbl 


265 MeV 


351 MeV 


Xf>2 


< 192 MeV 


213 MeV 



In Fig. [10] we show results for the real and imaginary parts of the binding energy 
of the T(ls) and the Y(2s)P^"^'lf the real part of the binding energy is positive, 
then the state is bound. If the real part of the binding energy is negative, then 
the state is unbound. The imaginary part of the binding energy gives information 
about the decay rate of the state in question. To see the exact relationship we can 
compute the quantum mechanical occupation number as a function of proper time 

n v{r) = (</)*(r,x)^(r,x)) , 

= (C(ro,x)^(r ,x))e 23 [ iJ ]^-^, 

= n° e 23 [ £ ](r-r o)) (32) 

where in the last line we have identified n£ = (0J(to, x)<^i,(to,x)) . In order to 
connect this to the decay rate, T, we note that T is defined empirically through 
n v (t) = n° exp(— r(r — r )) so that we can identify T — —2^s[E]. As a rough guide 
one can identify a dissociation temperature with the point at which the real and 
imaginary parts of the binding energy become comparable. At this point the state 
decays within one real-time oscillation of the quantum state. 

In practice, however, the precise value of the disassociation temperature is of 
little importance since even below this temperature the state under consideration is 
undergoing in-medium decays. In practice, one must compute the dynamical decay 
rate based on the local temperature and momentum-space anisotropy of the plasma 
and solve for the fraction of states remaining at a particular point in the plasma as 
a function of proper-time and then integrate/average over space-time. In order to 
carry this step out it is necessary to have a model of the spatio-temporal dynamics 
of the QGP. We will return this issue after the next subsection. 



3.3. Sequential Suppression and Feed Down 

As mentioned above, as rough guide one can determine a rough dissociation tem- 
perature for each state by finding the temperature at which the real and imaginary 
parts of the binding energies become equal. As can be seen from Fig. [10] (and Ta- 
ble [2]), the dissociation temperature determined in this way is higher for the T(ls) 
state than it is for the T(2s) state meaning that the T(ls) will survive to higher 
temperatures. In Table [2] we list estimates of disassociation temperatures for the 
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Fig. 11. Spectroscopic diagram for charmonium states. Decay channels shown were taken from 
Ref. 11321 The bottom row shows the spin, parity, and charge conjugation quantum numbers asso- 
ciated with the states above it. 

J/tp, Xci, T(ls), T(2s), T(3s), \bi, an d Xb2 states obtained in Refs. I129l II 30! In 
those references the authors determined the real and imaginary parts of the binding 
energies for each of the states listed. This is necessary since if, for example, one is 
interested in the suppression of the J/tJj state, it is necessary to calculate the direct 
suppression of all states that can potentially decay into the J/ip- 

The complications implied by feed down can be inferred from Fig. [TT] in which 
we have illustrated the spectrum of charmonium states along with the decay chan- 
nels of each stateP^I W e d no t show extremely rare decays or possible feed downs 
from bottom hadrons. As can be seen from this figure, when computing J/tp sup- 
pression, for example, it is necessary to have information about the suppression 
of the Vt'(2s), x c {i ! 2,3}(lp)j an d h c (lp) states. A similar complication arises in the 
bottomonium sector. In vacuum the percentage of produced, e.g., J/ip states which 
come from a particular decay channel is called the feed down fraction. Many experi- 
ments have attempted measurements of the feed down fractions for the charmonium 
and bottomonium systems. In Table [3] we present feed down fractions for J/tp and 
T production. The quoted numbers internal to the charmonium system are world 
averages presented in Ref. 11331 The J/ip bottom hadron (6-hadron) feed down frac- 
tion is estimated from Ref. 134. The direct J/ip production fraction is estimated 
by requiring that the sum of the listed feed down fractions be equal to one. The 
bottomonium feed down fractions are quoted from Ref. 11351 

As we can see from Table El in order to make quantitatively reliable calcula- 
tions of quarkonium suppression, it is necessary to calculate the suppression of the 
first few excited states. The direct suppression of each state (excited or ground) 
should then be multiplied by the corresponding feed down fraction and summed 
appropriately in order to obtain the inclusive suppression factor. 
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Table 3. Feed down fractions for (left) J/ip and (right) T production. Quoted numbers internal 
to the charmonium system are world averages presented in Ref. |133l The J/ip bottom hadron 
(6-hadron) feed down fraction is estimated from Ref. 134 The direct J/ip production fraction is 
estimated by requiring that the sum of the listed feed down fractions equal one. Bottomonium 
feed down fractions are quoted from Ref. 11351 



J/ip feed down fractions 


Mechanism 


% ± Error 


Direct Production 


41 ± 17 


*(2s) decay 


8.1 ± 0.3 133 


Xd decay 


25 ± 5.0 133 


6-hadron decay 


8.1 ± 3.2 134 



T(ls) feed down fractions 


Mechanism 


% ± Stat ± Sys 135 


Direct Production 


50.9 ± 8.2 ± 9.0 


T(2s) decay 


10.7 ± 7.7 ± 4.8 


T(3s) decay 


0.8 ± 0.6 ± 0.4 


Xbl decay 


27.1 ± 6.9 ± 4.4 


Xb2 decay 


10.5 ± 4.4 ± 1.4 



3.4. Dynamical Models for Quarkonium Production 

Phcnomcnological modeling of quarkonium suppression in the QGP is complicated 
by the fact that the matter generated in a heavy ion collision only exists in the 
plasma state for a few fm/c and during this very short time the properties of the 
QGP, such as the local temperature, rapidly evolve. At the LHC energies, for ex- 
ample, in central collisions the QGP is formed at r = To ~ 0.2 fm/c after the 
initial nuclear impact and has a temperature on the order of 600 MeV; however, 
the plasma rapidly cools to a freeze-out temperature T = Tf ~ 150 MeV within 
approximately 7 fm/c. For non-central collisions, the initial temperature decreases 
and, as a consequence, the lifetime of the plasma becomes shorter. In addition, the 
matter generated at all centralities has (on average) a spatial profile with the high- 
est temperatures in the center of the overlap region and lower temperatures towards 
the edges. In order to make quantitatively reliable predictions for quarkonium sup- 
pression, it is necessary to use dynamical models which can accurately describe the 
full spatio-temporal evolution of the QGP. 

Different groups have employed a variety of different models including basic 
"fireball" models (see e.g. Refs. 1X361 1137[) . ideal (2+l)-dimensional hydrodynamics 
including fluctuations (see e.g. Refs. 11381 11391 1140p . and (resummed) anisotropic 
viscous hydrodynamics EHHlUJIJCEII] At the very least, these dynamical models pro- 
vide information concerning the spatio-temporal evolution of the local QGP energy 
density; however, they can additionally provide information about, for example, 
the local flow velocity of the medium, local flow velocity of the heavy quark states 
themselves, the degree of momentum-space isotropy of the QGP, etc. This dynamical 
information must be integrated in proper-time in order arrive at the final predictions 
for the level of suppression for each quarkonium state under consideration. 

There are also other models that try to explain quarkonium production and 
suppression in heavy ion collisions. Most of these models ignore possible plasma 
anisotropies and assume complete local thermal equilibrium. Some of the models 
currently available are based on statistical recombinationf^2l statistical recombina- 
tion and dissociation rates^^ESIEI] (this is the so-called two component model); or 
sequential melting For loosely bound quarkonium states, as well as for unbound 
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QQ , a more realistic model has been developedP^S This incorporates a realistic 
modeling of the expanding hot medium and a realistic potential, and it is relevant, 
in particular, for the understanding of charmonium production in heavy ion colli- 
sions since, as we have seen in the previous section, no charmonium states can exist 
deep in the deconfined state. 

The bulk evolution of the matter produced in heavy-ion collisions is reasonably 
well-described by hydrodynamics (see e.g. Ref. 11481 for a recent review). The large 
heavy quark mass makes it possible to model its interaction with the medium by 
Langevin dynamics Such an approach successfully describes the anisotropic flow 
of charm quarks observed at RHKjSUllSQI ( se e also the review in Ref. 11511 and ref- 
erences therein). Potential models have shown that, in the absence of bound states, 
the QQ pairs are correlated in spaceP^E3 This correlation can be modeled classi- 
cally using Langevin dynamics, including a drag force and a random force between 
the Q (or Q) and the medium, as well as the forces between the Q and Q described 
by the potential. It was shown that a model combining an ideal hydrodynamic ex- 
pansion of the medium with a description of the correlated QQ pair dynamics by the 
Langevin equation can describe charmonium suppression at RHIC quite wellP^S I n 
particular, this model can explain why, despite the fact that a deconfined medium is 
created at RHIC, there is only a 40 — 50% suppression in the charmonium yield. The 
attractive potential and the finite lifetime of the system prevents the complete de- 
correlation of some of the QQ pairs Once the matter has cooled sufficiently, these 
residual correlations make it possible for the Q and Q to form a bound state. The 
model was further developed to include a realistic viscous hydrodynamic description 
of the medium 1152J Production of charmonium states resulting from recombination 
of originally uncorrelated cc pairs, so-called off-diagonal recombination, can be also 
calculated in this frameworkP^K^ The off-diagonal recombination is sub-dominant 
at RHICp^but is the dominant mechanism of J/tp production at the LHCp^In 
particular, the calculation in RefP^J W as able to explain the very mild centrality 
dependence of the J/tp production at LHC. 

3.5. Cold Nuclear Matter Effects 

The goal of studying quarkonium suppression in nucleus-nucleus (AB) collisions 
is to use them as a probe to determine information about the QGP; however, in 
practice one has to be careful since propagation in the QGP is not the only way in 
which quarkonium yields can be affected. As mentioned in the introduction, already 
in proton-nucleus (pA) collisions one sees modifications of quarkonium production. 
In order to compute the suppression of quarkonium states in a QGP one needs to 
know the initial spectrum of states including possible modifications from scatter- 
ing, interference, absorption, and energy loss while the precursor states propagate 
through the nucleus itself. Such modifications are lumped together into what are 
known as cold nuclear matter (CNM) effects. 

There are several CNM effects of importance. The first to consider is modification 
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of nuclear parton distribution functions (nPDFs). Due to this effect, in different 
kinematical regions quarkonium production can be suppressed (shadowing) and 
enhanced (anti-shadowing) relative to proton-proton (pp) collisions (we will lump 
both together as "shadowing" in the remainder). The key inputs necessary are the 
quark and gluon nPDFs. The quark nPDFs are fairly well-constrained by nuclear 
deep-inelastic scattering (nDIS) experiments; however, the gluon distributions must 
be inferred from scaling violations and conservation laws. 

In order to compute CNM shadowing, typically parton distribution functions 
such as CTEQpH are used in conjunction with different shadowing parameteriza- 
tions, e.g. EKS98P S * 15 3 nDSgP 3 HKN,™ EPSOS,^ or EPS09.^ Unfortunately, 
there is a sizable variation of the resulting estimates depending on the choice of shad- 
owing parameterization and in some cases even within the uncertainties quoted by 
an individual parameterization, e.g. EPS09. In a recent work, a study of CNM shad- 
owing effects on J/ijj and T production at the top LHC energy of ^/snn = 7 TeV 
was performedP^ The analysis found that for central Pb-Pb collisions at y = 0, the 
nuclear suppression factor due to CNM shadowing effects alone was on the order of 
Raa — 0.5 ±0.4 for the J/tp and Raa — 0.65 ±0.3 for the T, with the uncertainties 
estimated by varying the EPS09 parameters within their stated uncertainties. The 
central values indicate that there could indeed be sizable CNM shadowing effects on 
both the J/ip and T at LHC energies; however, more experimental and theoretical 
work is needed to further eliminate uncertainties in the nPDFs. In this context we 
also point to Ref. in which an analysis of CNM effects on T production was 
performed, and found at central rapidity shadowing on the order of 10% in p-Pb 
collisions at LHC energies. 

Another CNM effect which must be taken into consideration is the nuclear ab- 
sorption of quarkonium bound states after they form. In general the effective ab- 
sorption cross section depends on the collision energy and the rapidity of the ob- 
served state. At sufficiently high energies, one expects negligible absorption cross 
sections since the formation times for these states are higher than the amount of 
time they spend traversing the target. Early works assumed a constant absorption 
cross section which might be different for singlet and octet states EMilSSI however, 
more recent studie^SESl show that the effective absorption cross section decreases 
with increasing energy which seems to be in agreement with recent experimental 
findings! 171 ^ 172 ! In addition, we once again point to RefP^I which found that for the 
T the survival probability is close to unity, implying a small effective absorption 
cross section for bottomonium states. In this context, we mention that one also 
has to consider energy loss during the traversal of the target. For a discussion of 
different energy loss models and their impact on quarkonium production we refer 
the reader to Ref. 11681 

Finally, we mention that at very high energies the small- ir partonic wave function 
of a heavy ion is expected to be dominated by a highly-occupied saturated gluonic 
state called the color glass condensate (CGC)P^The CGC is expected to play an 
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important role in the initial state suppression of charmonium production, since the 
saturation scale Q s is on the same order of magnitude as the charm quark massP^ 
This effect could be particularly important for forward J /if) production!^! The 
effect is expected to be smaller for bottomonium states due to their larger mass; 
however, it should be taken into account nonetheless. The J/ip nuclear modification 
factor was studied in Refs. 11761 and 11771 in the CGC framework, and it has been 
shown that RdA can be well-described in this approachP^I 

4. Conclusions and Outlook 

In this paper we have attempted to make a compact review of (a) the theory 
of quarkonium states in the quark gluon plasma and (b) the phenomenology of 
quarkonium production in relativistic heavy ion collisions. There has been significant 
progress since the original predictions of Matsui and Satz. Thanks to first princi- 
ples calculations available from lattice QCD it is possible to make non-perturbativc 
studies of the behavior of the heavy quark potential and quarkonium spectral func- 
tions. Additionally, there has been significant progress made in recent years within 
the effective theory approach at zero and finite temperature allowing a firm ba- 
sis for systematically improvable potential model calculations. Using the effective 
theory approach, somewhat unanticipated imaginary-valued contributions to the 
quarkonium potential have been identified in recent years. These imaginary-valued 
contributions can be shown to be directly related to the in-medium decay of quarko- 
nium states. Depending on the QGP energy density and state under consideration, 
the implied in-medium widths are on the order of 20-100 MeV, which is orders of 
magnitude higher than the corresponding vacuum widths. 

Such large widths imply in-medium heavy-quarkonium lifetimes on the order of 
fm/c, thereby causing a blurring of the heretofore assumed pseudo-sharp thresh- 
olds for the dissolution of different quarkonium states, commonly referred to as 
dissociation temperatures. As a result, it is necessary to include the effect of such 
imaginary- valued contributions to the potential in phenomenological models and a 
few authors have already begun to implement the necessary dynamical models. In 
addition, phenomenological models of quarkonium suppression are seeing an across- 
the-board improvement in terms of the underlying dynamical models being used. 
Most practitioners are moving towards using some variant of hydrodynamics (ideal, 
viscous, or anisotropic) which allows the dynamics of collective flow and heavy 
quarkonium suppression to be described simultaneously. In addition, advances have 
been made in the modeling of the process of recombination, even including effects of 
QQ spatial correlations in a statistical approach utilizing Langevin dynamics. Need- 
less to say, further work is necessary but in recent years many theoretical groups 
have been contributing to key advances in this area. 

Looking forward, as the question of quarkonium suppression becomes more quan- 
titative, it will become increasingly important to further elucidate cold nuclear mat- 
ter effects on quarkonium production in heavy ion collisions. The pA runs at LHC 
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will provide much-needed input which will hopefully help to constrain models of cold 
nuclear matter effects. Additionally, further refinements of the space-time models 
of the quark gluon plasma are necessary including the incorporation of fluctuations 
in the viscous hydrodynamical treatments, further investigation of the impact of 
momentum-space anisotropics, effects of finite velocity of quarkonium states rel- 
ative to the medium, etc. Through the concerted effort of both theoretical and 
experimental collaborations, we hope to see continued progress towards a better 
understanding of the production and dynamics of quarkonium in the quark gluon 
plasma. 
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